Kmeans e mais exemplos

require(GGally, quietly = TRUE)
require(reshape2, quietly = TRUE)
require(tidyverse, quietly = TRUE, warn.conflicts = FALSE)

package ‘tidyverse’ was built under R version 3.3.2package ‘tibble’ was built under R version 3.3.2package ‘tidyr’ was built under R version 3.3.2

library(ggfortify)
library(cluster)
library(ggdendro)
library(broom)
theme_set(theme_bw())
source("github-lib.R")

dw <- load_github_wide()
summary(dw)

   repository_language   ForkEvent       IssuesEvent       PushEvent      
 ActionScript:  1      Min.   : 1.000   Min.   : 1.000   Min.   :  1.000  
 Ada         :  1      1st Qu.: 1.509   1st Qu.: 3.437   1st Qu.:  7.052  
 Agda        :  1      Median : 2.083   Median : 4.750   Median :  9.314  
 ANTLR       :  1      Mean   : 2.454   Mean   : 7.311   Mean   : 10.921  
 Apex        :  1      3rd Qu.: 2.913   3rd Qu.: 7.269   3rd Qu.: 10.602  
 AppleScript :  1      Max.   :18.000   Max.   :63.000   Max.   :154.250  
 (Other)     :121                                                         
   WatchEvent    
 Min.   : 1.000  
 1st Qu.: 2.000  
 Median : 3.007  
 Mean   : 3.725  
 3rd Qu.: 4.636  
 Max.   :13.471  
                 

dw %>% 
    select(-repository_language) %>% 
    ggpairs()

 plot: [1,1] [====---------------------------------------------------------]  6% est: 0s 
 plot: [1,2] [========-----------------------------------------------------] 12% est: 1s 
 plot: [1,3] [===========--------------------------------------------------] 19% est: 2s 
 plot: [1,4] [===============----------------------------------------------] 25% est: 1s 
 plot: [2,1] [===================------------------------------------------] 31% est: 1s 
 plot: [2,2] [=======================--------------------------------------] 38% est: 1s 
 plot: [2,3] [===========================----------------------------------] 44% est: 1s 
 plot: [2,4] [==============================-------------------------------] 50% est: 1s 
 plot: [3,1] [==================================---------------------------] 56% est: 1s 
 plot: [3,2] [======================================-----------------------] 62% est: 1s 
 plot: [3,3] [==========================================-------------------] 69% est: 1s 
 plot: [3,4] [==============================================---------------] 75% est: 0s 
 plot: [4,1] [==================================================-----------] 81% est: 0s 
 plot: [4,2] [=====================================================--------] 88% est: 0s 
 plot: [4,3] [=========================================================----] 94% est: 0s 
 plot: [4,4] [=============================================================]100% est: 0s 
                                                                                         

# XML e Bluespec têm mais de 50 pushes / repositório e 
# outras linguagens têm também números estranhos. Filtraremos.
dw <- dw %>% 
  filter(PushEvent < 50, IssuesEvent < 50, ForkEvent < 18)

As variáveis são bastante assimétricas e concentradas em pequenos valores. Transformá-las para log ajuda na visualização.

summary(dw2.scaled)

   repository_language   ForkEvent         IssuesEvent         PushEvent      
 ActionScript:  1      Min.   :-1.80630   Min.   :-2.22415   Min.   :-5.1947  
 Ada         :  1      1st Qu.:-0.81777   1st Qu.:-0.46527   1st Qu.:-0.4681  
 Agda        :  1      Median :-0.01687   Median :-0.02479   Median : 0.1830  
 ANTLR       :  1      Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.0000  
 Apex        :  1      3rd Qu.: 0.75644   3rd Qu.: 0.55363   3rd Qu.: 0.4992  
 AppleScript :  1      Max.   : 2.56203   Max.   : 2.76710   Max.   : 2.7606  
 (Other)     :115                                                             
   WatchEvent         cluster         
 Min.   :-1.88615   Length:121        
 1st Qu.:-0.72590   Class :character  
 Median :-0.04339   Mode  :character  
 Mean   : 0.00000                     
 3rd Qu.: 0.68503                     
 Max.   : 2.35217                     
                                      

dists = dw2.scaled %>% 
    column_to_rownames("repository_language") %>% 
    dist(method = "euclidean")
hc = hclust(dists, method = "ward.D")
plot(hc, cex = .6)

plot(hc, hang = -1)
n_clusters = 4
rect.hclust(hc, k=n_clusters)

dw2 <- dw2 %>% 
    mutate(cluster = hc %>% 
               cutree(k = n_clusters) %>% 
               as.character())
dw2.scaled <- dw2.scaled %>% 
    mutate(cluster = hc %>% 
               cutree(k = n_clusters) %>% 
               as.character())
dw2.long = melt(dw2.scaled, id.vars = c("repository_language", "cluster"))
hc %>% 
    cutree(k = n_clusters) %>% 
    silhouette(dists) %>% 
    plot(col = RColorBrewer::brewer.pal(n_clusters, "Set2"))

dw2.long %>% 
    ggplot(aes(x = variable, y = value, group = repository_language, colour = cluster)) + 
    geom_line(alpha = 0.4) + 
    facet_wrap(~ cluster) 

library(plotly)
p <- dw2 %>%
    plot_ly(type = 'parcoords',
            line = list(color = ~cluster),
            dimensions = list(
                #list(range = c(1, 4), label = "cluster", values = ~cluster),
                list(range = c(0, 4),
                     label = 'Sepal Width', values = ~ForkEvent),
                list(range = c(0, 4),
                     constraintrange = c(5,6),
                     label = 'Sepal Length', values = ~IssuesEvent),
                list(range = c(0, 4),
                     label = 'Petal Width', values = ~PushEvent),
                list(range = c(0, 4),
                     label = 'Petal Length', values = ~WatchEvent)
            )
    )
p

k-means

table(km %>% pull(cluster))

Error in function_list[[k]](value) : 
  não foi possível encontrar a função "pull"

Qual seria um bom valor de k? Uma medida comumente usada no kmeans é comparar a distância (quadrática) entre o centro dos clusters e o centro dos dados com a distância (quadrática) entre os pontos todos nos dados e o centro dos dados. Aqui o centro dos dados é um ponto imaginário na média de todas as variáveis. Calculamos a distância do centro de cada cluster para o centro dos dados e multiplicamos pelo número de pontos nesse cluster. Somando esse valor para todos os clusters, temos betweenss abaixo. Se esse valor for próximo do somatório total das distâncias dos pontos para o centro dos dados (totss), os pontos estão próximos do centro de seu cluster. Essa proporção pode ser usada para definir um bom valor de k. Quando ela para de crescer, para de valer à pena aumentar k.

set.seed(123)
explorando_k = tibble(k = 1:15) %>% 
    group_by(k) %>% 
    do(
        kmeans(select(dw2.scaled, -repository_language), 
               centers = .$k, 
               nstart = 20) %>% glance()
    )
explorando_k %>% 
    ggplot(aes(x = k, y = betweenss / totss)) + 
    geom_line() + 
    geom_point()

---

K-means

filmes = readr::read_csv("dados/filmes-scarlett-johanssson.csv")

Parsed with column specification:
cols(
  RATING = col_double(),
  TITLE = col_character(),
  CREDIT = col_character(),
  `BOX OFFICE` = col_double(),
  YEAR = col_integer()
)

Mais um exemplo

O dataset ruspini é clássico para ilustrar agrupamento.

Hierárquico

dists = dist(rs, method = "euclidean")
hc = hclust(dists, method = "ward.D")

plot(hc, hang = -1, cex = 0.8)

rect.hclust(hc, k=4)

rs$cluster = factor(cutree(hc, k=4))

ggplot(rs, aes(x = x, y = y, colour = cluster)) + 
  geom_point(size = 3) 

rs$cluster = factor(cutree(hc, k=8))
ggplot(rs, aes(x = x, y = y, colour = cluster, label = cluster)) + 
  geom_point(size = 2) + 
  geom_text(hjust = -.1, vjust = 1) + 
  xlim(0, 150)

plot(silhouette(cutree(hc, k = 4), dists))
plot(silhouette(cutree(hc, k = 6), dists))

#heatmap(as.matrix(dw2[,1:4]), Colv=F, scale='none')
#hc.data <- dendro_data(hc)
#ggdendrogram(hc.data, rotate = TRUE) + 
  #labs(title = "Agrupamento de Rustini")

km <- kmeans(rs, centers=4, nstart=10)
km

autoplot(km, data = rs)

autoplot(km, data = rs, frame = TRUE)